#include <iostream>
#include <stdio.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <vector>

#include "Timer/timer.01.h"
#include "Vector/vector.03.h"
//#include "Sorter/sorter.01.h"

using namespace std;

class ICC {
	public:
		
		static const double pi = 3.14159265358979323846;
		
		double m;
		
		typedef unsigned int uint;
		typedef long unsigned int luint;
		
		virtual double M(double r) {
			return 2*pi*r*r;
		}
		
		virtual double R(double m) {
			return sqrt(m/(2*pi));
		}
		
		
		template <typename T>
		vector<T> getr(uint N, double r1, double r2) {
			
			vector<T> r(N);
			
			double M1 = M(r1);
			double M2 = M(r2);
			
			double DM = M2 - M1;
			
			m = DM/N;
			
			printf("M2-M1= %g \n", DM);
			
			for (uint i=0; i<N; i++)
				r[i] = R( M1 + DM*(double)rand()/RAND_MAX );
			
			return r;
		}

		
		template <typename T>
		vector< Vector<T> > getR(uint N, double r1, double r2) {
			
			
			vector<T> r = getr<T>(N, r1, r2);
			
			vector< Vector<T> > R(N);

			for (uint i=0; i<N; i++) {
				
				T cosa1 = 1-2*(double)rand()/RAND_MAX;
			
				T a1 = acos(cosa1); 
				T a2 = 2*pi*(double)rand()/RAND_MAX; 
			
				R[i].x = r[i]*sin(a1)*cos(a2);
				R[i].y = r[i]*sin(a1)*sin(a2);
				R[i].z = r[i]*cosa1;
			}
			
			return R;
		}
		

		
		
		template <typename T>
		vector< Vector<T> > getHernquistR(uint N, double a, double r1=0, double r2=1e50) {
		
			class MyICC : public ICC {
			public:
				double a;
				double M(double r) {
					return r*r/pow(r+a, 2);
				}
		
				double R(double m) {
					return a /( 1/sqrt(m) - 1);
				}				
				
			} icc;
			
			icc.a = a;
			return icc.getR<T>(N, r1, r2);
			
		}
		
	
};
